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ABSTRACT 

Precessing accretion discs have long been suggested as explanations for the long periods 
observed in a variety of X-ray binaries, most notably HerX-l/HZ Her. We show that an 
instability of the disc's response to the radiation reaction force from the illumination 
by the central source can cause the disc to tilt out of the orbital plane and precess in 
something like the required manner. The rate of precession and disc tilt obtained for 
realistic values of system parameters compare favourably with the known body of data 
on X-ray binaries with long periods. We explore other possible types of behaviour than 
steadily precessing discs that might be observable in systems with somewhat different 
parameters. At high luminosities, the inner disc tilts through more than 90 degrees, 
i.e. it rotates counter to the usual direction, which may explain the torque reversals 
in systems such as 4U 1626—67. 
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1 INTRODUCTION 

Soon after the discovery of HerX-1 (Tananbaum et al. 1972) 
it was found that beside the 1.24-sec pulse period and the 
1.7-day orbital period (marked by eclipses of the neutron 
star by the companion) the system displayed another ('long' 
or 'third') period: a 35-day cycle of X-ray on and off states. 
Photometric variations showed a period of 1.6 days, the beat 
period between the orbital and long periods. A phenomeno- 
logical 'clockwork' model in which the system contains a 
retrogradely precessing disc was found to explain the pho- 
tometry quite well. The photometric variations are due to 
the varying amount of area of the disc facing us, as well as 
variations in the fraction of the star that is eclipsed when the 
disc passes in front of it and the variable X-ray illumination 
of the star (Gerend & Boynton 1976, Deeter et al. 1976). The 
behaviour is not strictly periodic: there are 'anomalous' on- 
and off-states in X rays, and the period of the cycle varies 
by 5-10% over time (Ogelman et al. 1985). With time, two 
more systems were found with nearly strictly periodic long 
variations: SS 433 has jets which precess around the sky ev- 
ery 164 days, and therefore presumably so does the under- 
lying disc, and the direction of precession is retrograde (see 
Margon 1984). LMCX-4 has a 30.4-day long period, which 
resembles that of HerX-1 in many ways, except that the 
amplitude is much lower due to the fact that the companion 
star is much brighter (Heemskerk et al. 1994). 



The success of the clockwork model quickly led to the 
general acceptance of the basic hypothesis that the long pe- 
riods are due to precessing tilted accretion discs (Katz 1973, 
Roberts 1974). It appears easy to make a disc in a binary 
precess, since the companion star exerts forces on the disc 
which, when averaged over the orbit, lead to a constant rate 
of retrograde precession of a tilted annulus in the disc. The 
problem is to understand how the disc as a whole can pre- 
cess at one rate, since the precession frequency is a strong 
function of radius in the disc, and how it remains tilted, i.e. 
avoids sinking back into the orbital plane on a viscous time 
scale. Roberts proposed that the disc is slaved to a compan- 
ion whose spin is misaligned with the orbit (Roberts 1974). 
But that is difficult on general grounds: all neutron star bi- 
naries have eccentric orbits soon after the neutron star has 
formed, due to the effects of the supernova explosion. But 
now most (notably HerX-1) have orbits that are accurately 
circular, indicating that tidal forces have been strong enough 
to damp this initial eccentricity. It is hard to see how they 
would at the same time not have largely damped out a spin- 
orbit misalignment (Petterson 1977), since the spin angular 
momentum is much less than the orbital angular momentum 
for such binary systems. 

Beside the three very regular and clear third periods 
mentioned above, there are many systems where a less well- 
defined variability exists. Some are quite regular but of more 
uncertain nature, e.g. CygX-1 (294 d; Priedhorsky, Terrell 
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& Holt 1983) and 4U 1820-30 (176 d; Priedhorsky & Ter- 
rell 1984a). Others are simply less regular, such as CenX-3 
(120-160 d; Priedhorsky & Terrell 1983). In this latter cate- 
gory numerous objects have recently been added with well- 
sampled light curves obtained with the Rossi X-ray Tim- 
ing Explorer (RXTE). CygX-2 (78 d; Wijnands, Kuulkers & 
Smale 1996) and X 2127+119 (37 d; Corbet, Peele & Smith 
1997) have periods that seem fairly stable, but the ampli- 
tude and light curve shape are variable so the periodicity 
cannot always be seen equally well. SMCX-1, on the other 
hand, has a fairly stable amplitude but a period that slowly 
decreases from over 60 d to under 50 d in the first 600 days 
of RXTE data (Wojdowski et al. 1997). 

Moreover, there appears to be some evidence of discs 
precessing progradely. First, there has been a detection of 
an 11.2-day period in the light curve of CygX-2 (Holt et al. 
1976). This is the beat period between the 78-day long pe- 
riod and the 8.9-day orbital period, but only if the disc pre- 
cesses progradely. Note that this cannot be due to tidally 
forced precession. In X 1916—053, the optical period is 0.9% 
longer than the X-ray period, and if the beat period of 3.8 
days between the two is interpreted as the period of disc pre- 
cession, it would most easily be for a progradely precessing 
disc. 

It was noted by Petterson (1977) and Iping & Petterson 
(1990) that a sufficiently strong illumination from the centre 
could cause a disc to maintain a warped, tilted shape. While 
probably qualitatively correct, their results suffer from the 
use of an incorrect equation of disc evolution (see Papaloizou 
& Pringle 1983). An accretion disc is indeed unstable to 
tilting and warping due to radiation reaction forces when 
the luminosity of its central source exceeds a critical value 
(Pringle 1996, 1997). We examine here the consequences of 
this instability for the behaviour of discs in X-ray binaries. 
We show that the instability provides a mechanism for sus- 
taining a tilted disc and making it precess with a period 
that agrees well with the observed long periods in X-ray bi- 
naries such as the famous 35-day period of HerX-1. They 
also provide the possibility of both prograde and retrograde 
precession. There are also non-steadily precessing solutions 
with time- varying tilt. We explore these in Section ^| after 
discussing the numerical solution method to the basic equa- 
tion, which is derived in Section ^. Then we apply the results 
to some of the known long periods in X-ray binaries and to 
some similar systems that do not have observed long peri- 
ods (Section k|). Finally, we discuss some implications and 
limitations of our findings (Section ^ and summarise our 
conclusions (Section^. 



2 EQUATION OF MOTION FOR AN 
IRRADIATED DISC IN A BINARY 

The accretion disc is assumed to be thin and Keplerian and 
is divided up into annuli that interact with each other via 
viscous forces, which has the advantage that only the evo- 
lution on the longer viscous time scale needs to be followed. 
The relevant equation of motion reads 
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The first three terms on the righthand side were derived by 
Papaloizou & Pringle (1983, eq. 2.4) and further discussed 
by Pringle (1992, eq. 2.8). They conserve angular momentum 
exactly. The independent variables are {R,t) and we work 
in a non-rotating frame centred on the accreting point mass 
M (with the z axis normal to the orbital plane), t is the 
unit vector (cos 7 sin (3, sin 7 sin (3, cos /3) normal to the disc 
at radius R, so /3(R, t) is the local tilt angle of the disc plane 
and *y(R,t) the azimuth of the tilt. E is the disc surface 
density, and L is the angular momentum per unit area. The 
assumption of Keplerian rotation implies L = E^V GMR. 
v\ is the viscosity associated with the (R, (p) shear, whereas 
the viscosity associated with the (R, z) shear that damps the 
misalignment between neighbouring annuli is The term 
involving G is the effect of the torque induced by irradiation 
from the central object: 
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where the dimensionless vector is an integral along the 
annulus that describes the purely geometric part of the 
torque (see Pringle 1996, 1997). £* is the luminosity of the 
central point source, which is assumed to radiate isotropi- 
cally, and each element of the disc that is illuminated by the 
central source is assumed to re-radiate the incident radiation 
isotropically, causing a pressure perpendicular to the local 
disc surface. 

The last term in Eq. (|l| is the forced external precession 
due to the companion star tide. Let a companion of mass 
M c be located in the XY plane at a = (a cos ip, a sin ip, 0). 
To lowest order in R/a the torque on an annulus of width 
dR can be obtained by integrating the force along an unper- 
turbed circular orbit. The force on an element of the annulus 
at position x is 



dF c 



GM c (a - x) 



ERdRd<p. 



(3) 



The torque on the annulus is then simply obtained by in- 
tegrating the elemental torque contributions along it. Here 
we expand the force in powers of R/a and only retain the 
lowest non- vanishing order of the torque: 



xx dF c = — GM c Edi?^- 

2 a J 



sin 2/3 sin ip cos(7 — ip) 
—sin 2/3 cos ip cos(7 — ip) 
sin 2 /3 sin 2(7 — ip) 



(4) 



Note that the only force taken into account is the gravity 
of the companion. Since our frame is centred on M (but 
not rotating) and thus revolves around the centre of mass 
of the binary, there is a centrifugal force as well. It is not 
hard to show that it contributes no net torque on an annulus. 
Comparison with the expression for L leads to the expression 
for the precession frequency: 
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sin/3cos(7 — ip) 



In this paper, we shall make the approximation that the 
orbital period of the binary is small enough compared to the 
precession period that we may average the precession over 
the binary orbital period (i.e. over the angle if)). Averaging 
and cross-product are not interchangeable, of course, so we 
have to average the torque (Eq. (Q)) and then find a new 
expression for S~l p . The result is: 



f2p V = f2 p ep V = Sl p cos /3e z 



(6) 



(See also Katz et al 1982, Papaloizou & Pringle 1983). Note 
that the factor cos f3 is not present in studies that assume 
the inclination to be small. 

For numerical studies we cast the equation in dimen- 
sionless form. The dimensionless variables are r — R/Ro, 
a = E/Eo and r = t/to, with Ro, to and Eo values to 
be determined later. Then we can also set C — L/Lo = 
L/T,oVGMRo, v\ = nivio, and r\ = vi/v\. Then if we mea- 
sure time in units of the viscous time scale in the disc, i.e. 
to — Rq/v\o, we obtain: 
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where ui p = Sl p to, and the dimensionless strength F* of the 
radiation field is given by 

= 



L* Rp 
67rcEo-Rofio via 



(8) 



2.1 Analytic considerations 



The stability to warping of a disc described by Eq. (gj) can 
be analysed by linearising the disc tilt evolution (eq. 3.5 
of Pringle 1996) in /3. Then £ — (/3 cos 7, (5 sin 7, 1) and we 
can write its evolution as an equation for W = /3e 17 . We 
further assume that we are far from the disc edges, so that 
M oc j/iE = constant and Vr = viSl' /SI. Then with E oc 
^-f 1 oc v^ 1 we can write the following equation for W: 



dW 

dt 
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1 dR 2 + \ AR 
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where T = /12ttc'SR 2 S1. The local stability then follows 
by looking for solutions of the form W = Wo exp i(at + kR), 
with kR 2> 1. Substituting in Eq. <J^) this implies 

a = i[-Yk+\v 2 k 2 ] + Sl p + ^. (10) 

The complex part of a is identical to that of Pringle, so 
neither the addition of forced precession nor retention of 
the gradient terms affects the local linear stability of the 
disc. Thus, as before, growing modes occur only for Y > 
and < k < 2T /v-i and their line of nodes follows a leading 
spiral. At finite amplitudes the effect of forced precession is 
strong because Si p oc 7? 3//2 , so it will make the disc precess 



more rapidly towards the outside and thus tends to destroy 
the shape of the growing mode. 



3 NUMERICAL EXPERIMENTS 
3.1 Numerical Method 

The numerical method employed is that described in Pringle 
(1997), with changes made in order to model accretion discs 
in binary systems rather than in AGN. The most important 
changes which need to be made are with regard to tidal 
truncation of the disc at the outer edge, with regard to the 
input of mass at some radius other than the outer edge, and 
with regard to the tidally forced precession of the outer disc 
elements. 

We work initially with 40 grid points between r ln = 1 
and Tout = 40, and use an equally spaced rather than a 
logarithmic grid (Section 3.2). This was chosen so that the 
outer edge of the grid would be properly resolved, but has 
the disadvantage that the inner disc regions are not well 
resolved. However, because of this, the computation time 
required for each run is relatively short, and thus we were 
able to explore a lar ge volume of parameter space. We report 
below (Section 3.3) on a few runs carried out with a more 



extensive grid. The number of azimuthal zones was normally 
68. Increasing this value to 120, which we did in a few key 
runs, had no effect on the results. 

We achieve a zero torque inner boundary condition in 
the manner described in Pringle (1997) by removing mass 
over three zones centred on r — 3. We add material at a con- 
stant rate over 3 zones centred on radius r a( jd. The material 
is added to the disc with specific angular momentum vec- 
tor corresponding to the material being added in the orbital 
plane, and magnitude appropriate for the radius r a( jd. 

In order to take account of the tidal truncation of the 
disc at radius rude (< r out ), we need to set up a bound- 
ary condition which removes angular momentum at that ra- 
dius but does not remove mass. Since the main dependent 
variable employed is the angular momentum density C, this 
must of necessity employ some approximation. We find that 
the simplest method of achieving the desired effect is to cal- 
culate at each time step the amount of mass external to 
rtidej to set the contents of those zones to zero, and to add 
the mass to the zone just internal to r t id c , with the same 
specific angular momentum as the matter already in that 
zone. In this way it was intended to minimise the chance 
of the outer tidal truncation resulting in any (unphysical) 
instability. Indeed no such instabilities were observed. This 
boundary condition should correspond to the combined con- 
ditions of vr — and dl/dR = at r — r t id c - 

The forced precession caused by the tidal potential field 
is handled as follows. Each annulus is precessed at each time 
step about the vector (0,0,1) at a rate uj p — ui p or 3 ^ 2 . Thus 
we have ignored the factor of cos f3 present in the forced pre- 
cession rate (Equation 1.6). This could be simply included, 
but should be negligible for the kind of solutions we are look- 
ing for. The time step is limited to ensure that any annulus 
is not precessed by more than 27r/10 in any one time step. 
In practice the time step considerations for stability of the 
numerical method ensure that the precession rate is not the 
dominant factor in limiting the time step. 
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We assume a viscosity of the form: 

vi = v w r 3/4 , (11) 

which is chosen to approximate the viscosity dependence 
in such discs (see eq. Il5|) . We take Tj = 1. The disc is set 
up initially with a surface density appropriate to a steady 
disc with this viscosity and with the appropriate boundary 
conditions (viz. E = at r = 1, vr = at r = rude, and 
a dimensionless mass input rate of unity at r — r a dd). The 
disc is given an initial warp of the form of a prograde spiral 
which is such that j3 increases from zero at r = 1 to /? = 0.1 
at r = 40, with 7 increasing by an angle of one radian over 
the same distance. 

3.2 Numerical Results 

We have investigated the non-linear behaviour of irradiated 
discs in binary systems, taking shadowing fully into account. 
Our aim has been to discover what regions of parameter 
space (if any) give rise to discs with a finite, and (reasonably) 
steady tilt angle of the kind which might be relevant to ex- 
planations of systems like Her X-l/HZ Her. Thus our initial 
aim was to search for the kind of 'mode-like' behaviour dis- 
cussed by Maloney & Begelman (1997). We found, however, 
that the time-dependence of binary accretion discs gives rise 
to a richer and more varied behaviour than can be described 
by simple 'modes'. 
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Figure 1. The behaviour of the disc inclination at r = 26, just 
inside rt iHo, for different values of and using the small grid 
(Section p.2|). With increasing F+, the disc first is stable, then 
grows to a finite and constant tilt that is greater for higher values 
of F*. The oscillations at still higher F* are an artefact of the 
poorly resolved inner disc. 



8.2.1 r add = 10, r- ti do = 30, lj p0 = 

With an eye on the parameters of the binary systems we are 
interested in as shown in Table ^ we choose r a dd/r t idc = 
0.33. We initially set the forced precession rate to zero, and 
investigate what values of the disc luminosity (F+) give rise 
to relevant behaviour. In Figure |l| we show the behaviour 
of the inclination of the outer disc edge as a function of 
dimensionless time for various values of F*. We find that for 
F+ < 0.04 the disc irradiation is not strong enough to lead 
to instability and the disc flattens into the orbital plane. 

For F* — 0.045 the disc tilt settles down to a steady 
functional form, with /3 small at the inside, and /3' positive 
at most radii, reaching a value of j3 = 0.275 (corresponding 
to an angle of 16°) at the outside (i.e. at r = /"tide)- The 
disc is warped in the shape of a prograde spiral, as is to 
be expected for a disc unstable to self-irradiation (Pringle 
1996). At the same time the whole disc precesses steadily in 
a retrograde direction with a dimensionless period of P p = 
58. The direction of precession is retrograde because of the 
shape of the disc. The sign of determines which face of 
the disc is illuminated by the central source. A positive /3' 
leads to retrograde precession (Pringle 1996). We note that 
this contrasts with the AGN discs discussed in Pringle 1997. 
There the outer boundary condition ensured that the outer 
edge of the disc remained in the initial disc plane. Then if 
the inner disc is tilted, much of the disc is likely to have 
negative /3', and therefore precesses in a prograde fashion. 
In the binary case the edge of the disc is free in the sense 
that no torque acts there to change the tilt of the disc. 

Similar behaviour occurs when F* = 0.05, although in 
this case (3 at the outside settles down to the larger value of 
(3 = 0.44 (corresponding to an angle of 25°) and to a more 
rapid retrograde precession rate with a period P p = 48. 



When the effect of radiation is further increased, the 
disc behaviour is no longer steady. The behaviour of (5 at 
the outer radius is shown in Figure |l|. The behaviour is in 
the form of a limit cycle. We note here that for simplicity 
we have assumed that F t is constant throughout each run. 
However, in reality, since the radiation illuminating the disc 
comes directly from the accretion rate at the disc centre, F* 
is likely to vary with time when the disc displays such non- 
steady behaviour. Thus the actual disc behaviour is likely 
to be yet more complicated than we find here. 



3.2.2 r add = 10, r ti d c = 30, F* = 0.05 

In this section we describe the effect on the steady precessing 
solutions of adding retrograde forced precession of the form 
induced by tidal torques from a companion. The magnitude 
of the forced precession is described by the parameter ui p o 
(eq. ^). With zero forced precession, uj p o = 0, we have seen 
(Section |3.2.1 ) that the disc settles down in the shape of a 
prograde spiral, with j3 an increasing function of r, reaching 
(3 = 0.44 at the outer edge, and precessing steadily in a 
retrograde direction with period P p = 48. 

When a small amount of retrograde forced precession 
is added, the disc settles down to a similar behaviour as 
for lu p o — 0, but with the final value of /3(r t id e ) slightly 
increased, and the precession rate also increased. Thus we 
find that for w p0 = -0.0002, /3(r tidc ) = 0.54 and P p = 40; 
for lu p o — —0.0005, f3(r tidc ) = 0.63 and P p = 35, and for 
w p o = —0.001, /3 (rtidc) = 0.56 and P p = 29. However, as the 
forced precession rate is increased further the disc instability 
is removed. Thus, for u> p o = —0.002, we find that /3(r t id e ) 
tends to zero, and the disc settles down into the orbital 
plane. We suggest that this behaviour comes about because 
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in order to be unstable the disc must take the form of a 
progradc spiral (Pringle 1996). However strong retrograde 
forced precession, which is differential in the sense that the 
precession rate increases with radius (here oc r 3 ^ 2 ), tends 
to unwind the prograde spiral and thus acts to prevent the 
instability from occurring. 



3.2.3 



r-add 



10, rtide = 30, F* = 0.09 



Here we investigate the effect of adding forced retro grade 
precession on the solutions discovered in Section 



3.2.1 



which 

take the form of limit cycles. We therefore investigated so- 
lutions with F* = 0.09, and with uj p0 = -0.001, -0.002, 
—0.003 and —0.004. Increasing the forced precession rate 
seems to have little effect on the amplitude of the limit cy- 
cle, until at the value of —0.004 the instability is quenched 
altogether, and the disc settles into the orbital plane. How- 
ever the period of the limit cycle is affected, and decreases 
with increasing uj p o, from a value of 200 for lu p q = to a 
value of about 120 for ui p q — —0.003. While these limit-cycle 
solutions are partly an artefact of the small grid, the damp- 
ing of even these high-luminosity solutions illustrates the 
how powerful forced precession is in suppressing disc tilts. 



3.2.4 r- add = 20, r tidc = 30, uj pQ = 

One effect of adding material at r a dd = 10 in the previous 
sections was to help pin the disc towards the orbital plane 
at that radius. This helped to control the behaviour of the 
disc within that radius (see also Section 3.3), and to allow 
the outer regions of the disc (a factor of three in radius) to 
evolve freely with regard to tilt. Since the self-irradiation 
warping instability acts more strongly at larger radii, it was 
always the outer regions of the disc which responded most 
to the radiation flux from the central object. In this section 
we describe the effect of changing the radius at which mass 
is added to the disc from r a dd = 10 to r a dd = 20. This has 
two major effects. First, adding matter closer to the outside 
has the effect of tending to pin the disc into the orbital 
plane at the outside, whil e allowing the inner disc regions to 
tilt (see also Section 3.3). Second, a steady disc inside the 
radius r a dd has z/E constant at radii inside r a dd, whereas 
a steady disc with a vr = outer boundary condition has 
j/S oc r~ 1//2 for radii outside r at jd- Thus for a given accretion 
rate, since v is typically an increasing function of both £ and 
r, the outer parts of the disc become less massive as r a dd is 
decreased. Thus by increasing r a( jd from 10 to 20, we expect 
the radiation instability to require larger values of F*. 

For values of F* < 0.15, we find that the disc is stable 
to self-irradiation. For values of F* > 0.25 the inner parts 
of the disc turn over completely in the manner discussed for 
the AGN discs in Pringle 1997. For F* = 0.2 the disc does 
tend to a fairly steady configuration. However the shape of 
the disc now differs from those discussed above in the sense 
that the disc tilt f3 is largest at the inside, and decreases with 
radius (i.e. j3' is predominantly negative; Fig. ^). Moreover 
the inner and outer parts of the disc behave in quite different 
and independent manners. The inclination at all radii in the 
disc oscillates more or less in phase with a period of about 
Fine = 10. The innermost radii oscillate between j3 = 0.53 
and 0.57, and the outermost radii oscillate between j3 = 0.06 




radius 

Figure 2. Comparison of the variation of the disc tilt /3 with 
radius for mass addition well inside the outer radius (r add = 10, 
solid curve) and mass addition nearer the outer edge (r add = 20, 
dashed curves). The former, having /3' > in most of the disc, 
precesses retrogradely, whereas the latter precesses progradely. 
Note that when r add = 20 the disc is not a fixed-shape object 
precessing at a single rate (see text). The two curves span ap- 
proximately the range of shapes it displays. 



and 0.12. At radius r = 24, the oscillation in f3 varies be- 
tween and 0.14, and it is this radius which appears to sep- 
arate the inner and outer parts of the disc. The inner part of 
the disc (roughly those radii r < 24) has an inclination which 
decreases with radius, and which precesses in a prograde di- 
rection with a period of about F n = 12. The outer parts of 
the disc (roughly radii in the range 24 < r < 30) have (on a 
time average) an inclination /3 which increases with radius 
and precess in a retrograde direction with a period of about 
F out = 50. The precession in the inner regions proceeds at 
a more or less steady rate. However the azimuth of the tilt 
of the outer edge of the disc stays almost constant for each 
oscillation period of the inclination, and then jumps by a 
certain amount (which then gives the mean outer precession 
period). It should be noted that the oscillation frequency 
of the inclination is the sum of the moduli of the preces- 
sion frequencies of the inner and outer parts of the disc (i.e. 
F" 1 = Pr 1 + p- 



- 1 ) 

out / 



3.2.5 r add = 20, r tid c = 30, F* = 0.2 

We now investigate the effect of adding retrograde forced 
precession to the F* = 0.2 disc which was described in 
the previous section. As the magnitude of w p o increases 
from zero, the behaviour of the disc stays initially the 
same, except that the inclination decreases, the inner (pro- 
grade) precession period, Fi n , increases, and the outer pre- 
cession period, F out , decreases. Thus, in comparison with 
P(rin) = 0.55, F n = 12 and F out = 50 for lo p o = 0, we find 
that for w p0 = -0.0015, (3(r in ) = 0.46, F in = 14, F out = 20; 
for cu p0 = -0.004, f3(r in ) = 0.32, F in = 20, F out = 11; for 
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w p0 = -0.006, /3(r in ) = 0.2, P in = 26, P out = 6; and for 
w p0 = -0.007, /3(r ln ) = 0.16, P in = 40 and P out = 4.5. In 
addition, as the size of ui-po is increased, the radius which sep- 
arates the inner disc and outer disc behaviours (the radius 
at which the f3 oscillation passes through zero) moves out- 
wards until at ui p o = —0.006 it is at the outer edge to within 
the resolution of the grid. However when uj p o is increased 
further to —0.008 the whole disc now settles to a constant j3 
profile, which is positive at all radii, and has /3(n n ) = 0.12, 
and f3 decreasing with radius. The disc as a whole precesses 
in a prograde direction with period P p = 22. The shape of 
the disc is such that it has a prograde spiral inside r a dd, and 
a retrograde spiral outside r a dd- 

3.3 Numerical results on a more extensive grid 

Since all the numerical results above were computed on a 
grid which had limited resolution in the inner regions, we 
felt that it was necessary to explore the limitations of such 
a procedure. We use a logarithmic grid consisting of 80 grid 
points between Pi n = 0.136 and Ptidc = 30. We still add 
material over 3 grid points centred on P a dd = 10, and we 
remove material over 3 grid points centred on R — 0.155. 
In this manner the outer part of the grid between P a dd and 
Ptido is reasonably well resolved (15 grid points), and the 
inner disc region extends inward of P a dd by almost two or- 
ders of magnitude, so that the behaviour of the inner parts 
of the grid can now be modelled more accurately. 

Because the grid is now more extensive, and especially 
because the grid now extends to smaller radii (where vis- 
cous time scales are shorter) the computational run times 
are now longer by about two orders of magnitude. Thus we 
have limited ourselves to a few representative examples for 
comparison with the results in the previous section. We find 
that the instability sets in at values of P* which are larger 
by about a factor of two. This comes about because in a 
steady accretion disc with an inner boundary condition cor- 
responding to vanishing surface density (i.e. zero torque), 
it is the quantity i/E[l — (r/nn) 1 ' 2 ] which is constant with 
radius, rather than just i/E. Thus the effect of moving the 
inner boundary inwards is to increase the surface density of 
the outer disc by about a factor of two. 

Since the tidally induced precession does not seem to 
play a strong role for those discs which are relevant to the X- 
ray binaries we are interested in here, we have taken uj p — 0. 
In Figure ^ we show the behaviour of the disc inclination at 
r = 26 for P* = 0.09, 0.12, 0.135, 0.15, 0.175, 0.2 and 0.3. 
As can be seen the disc with P* = 0.09 is stable. The disc 
with P* = 0.12 eventually precesses in a steady fashion in 
a retrograde direction with a period of about 35, and the 
disc inclination at the outer edge settles down to a value of 
0.15 (i.e. about 8.5 degrees). The disc with P* = 0.15 settles 
down to a solution in which the inclination and oscillates 
with a period of about 15 and semi-amplitude 5 per cent, 
about a steady value for the inclination of 0.25 (14 degrees). 
The precession period is 30 for the outer disc. The inner disc 
has inclination near zero, with small oscillations of the same 
period as the outer disc (Figure ^) . The period of these os- 
cillations is the beat period between the inner and outer disc 
periods, as with the small grid. Their cause is simply that 
the region where the outer and inner disc join (just outside 
r a dd) tries to adjust to the tilts of both sides. This it can- 



not do simultaneously, of course, and the tilt of this zone 
oscillates between nearly zero and a finite value, depending 
on the relative phase between the inner and outer disc so- 
lutions. This oscillation is then communicated throughout 
the disc, with an amplitude that decreases away from the 
contact zone. 

For greater P*, the inner disc suddenly gets a substan- 
tial tilt as well, presumably because it too is now unstable 
(the sudden transition is caused by the fact that the in- 
ner radius of the unstable region decreases as the inverse 
square of P*, see Pringle 1997). The discs with P* = 0.2 
and 0.3 vary chaotically, with the tilt of the inner disc usu- 
ally greater than 90 degrees, and that of the outer disc less. 
The outer disc still precesses retrogradely on average, with 
roughly the same period as before (25-35), but with irreg- 
ularities superimposed. The azimuth of the inner disc tilt 
wanders irregularly without any long-term trend. 

In conclusion we find that the results obtained with 
the larger grid are fully consistent with the behaviour found 
with the cruder grid. Thus the results presented in Section 
3.2 are expected to be reasonably representative. The main 
differences appear to be that the larger grid has a larger 
surface density (for the reasons explained above) and so goes 
unstable for somewhat larger values of P*, that the range 
of P* for which the disc displays steady, or nearly steady, 
precessing tilted behaviour is somewhat larger, and that the 
tilt angles reached by the outer disc are somewhat smaller. In 
order to check the grid-dependence of these large simulations 
we further doubled the number of radial and azimuthal grid 
points in a few key simulations; no significant changes in the 
results were observed. 



3.4 Precession periods and stability 

In order to see how well simple estimates work, we now 
compare the precession periods obtained in the simulations 
with expected values. From Eq. (Q), one obtains a precession 
time scale estimate, taking g$/2n to be unity: 



Cr 

F,' 



(12) 



where r = R/Rq. Now C = ar 1 ^ 2 , where a = E/Eo, and 
thus at the outside edge of the disc, 



P p = 27rf rp = ^ ar 3/2 



(13) 



We now have to take the quantities from simulations to com- 
pare the left an d righ thand sides numerically. Using the re- 



sults from Sect. 3.2.1, we have, e.g. a numerical period of 
48 when P» = 0.05. At the outer edge of this simulation, 
r = 30.25 and a = 2.26 x 10~ 3 , giving a predicted period 
from Eq. ( |l3| ) of P p = 47, in very good agreement with the 
value seen in the simulation. For the run with P» = 0.045 
we have a = 2.76 x 10 -3 , and hence compute P p = 64, com- 
pared with the actual simulation value of 58. It is clear that 
the radiative time scale at the outer edge of the grid predicts 
the numerical precession period quite accurately. The high 
precision is doubtless somewhat fortuitous, since shadowing 
will tend to lengthen the period at the grid edge somewhat 
< 1), whereas the fact that the pattern speed must 
be an average over a finite range of radii in the disc tends to 
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Figure 3. The behaviour of the disc inclination at r = 26, just inside r t id c , for different values of F*. With increasing F*, the disc first 
is stable, then grows to a finite and constant tilt that is greater for higher values of F*. The oscillations around the stable level are due 
to the fact that the inner and outer disc precess in opposite directions once the unstable region in the disc includes enough of the disc 
inside r at jd, so that the solution is no longer stationary. At even greater F* the disc tilts through 90 degrees and the behaviour becomes 
chaotic. This is the regime previously discussed by Pringlc (1997) for AGN. 
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Figure 4. The behaviour of the disc inclination at r = 2.5, well inside r a( jd, for different values of F*. Note the larger tilts in comparison 
with the outer disc and the fact that at large F* the inner disc is mostly counter-rotating. 
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shorten it, because the radiative time scale becomes shorter 
at smaller radii. 

The critical number for stability is the ratio of radiative 
to viscous time scales at the outer disc edge. In the scaled 
equations, we have 



7 = *L = — 



^20 1/4 



(14) 



In analogy with Pringle (1997, eq.3.9) we can then define 
a critical value 7c rit, such that discs with 7 < 7 cr i t will be 
unstable to radiative warping. Since its precise value will 
depend on boundary conditions, it has to be evaluated from 
appropriate numerical simulations. There was a difference 
in the critical value of F t between the large and small grids 
we explored, and we surmised that this was due to their dif- 
ferent outer-edge surface densities. The evaluation of 7cr i t 
should be independent of this. For the small grids, the crit- 
ical value of F, is 0.045, and with a = 2.76 x 10~ 3 at 
r — 30.25, we deduce 7 cl -i t = 0.106. For the large grids, 
the smallest unstable F, is 0.09, and for this simulation 



a = 5.28 x 10 3 at r = 29.86, resulting in 7c; 



0.102. The 



excellent agreement between the two cases confirms that we 
understand the difference between the stability in the large 
and small-grid simulations. The fact that we find a three 
times smaller 7cr i t than Pringle (1997) did for simulations 
geared to AGN (i.e. without a tidal truncation on the disc) 
does demonstrate the need for re-determining 7 cr i t for dif- 
ferent types of system if one wants to go beyond order of 
magnitude estimates (see also Maloney, Begelman & Novak, 
1998). 



3.5 Summary of numerical results 

Since we have presented the results of a variety of numerical 
computations, using two different numerical grids, it would 
seem prudent at this stage to provide a brief summary of 
the basic results we feel that we have learnt from them. 

The results from both grids agree that for values of 
F* that are too small the disc is stable against the tilting 
instability. They further agree that for values of F* above 
the critical value for instability by not more than about 30 
per cent, the disc takes on a fixed shape, with inclination 
increasing to (at most) around 20 degrees at the outside 
(dependent on F+) and in the form of a prograde spiral, and 
precesses steadily in a retrograde direction (at a rate also 
dependent on F±). The extent of the outer part of the disc 
which is able to display this behaviour depends to some ex- 
tend on the radius at which mass is added to the disc. For 
values of F* slightly larger than this, the disc displays sim- 
ilar behaviour in the mean (i.e fixed disc shape and steady 
retrograde precession), but has an oscillatory behaviour su- 
perimposed upon this due to the fact that the inner and 
outer disc precess in opposite directions, which means that 
in the zone where they meet there is no stationary solu- 
tion and the tilt oscillates at the beat period between the 
inner and outer disc precession periods. This oscillation is 
communicated throughout the disc, with an amplitude that 
decreases away from the contact region. 

We have used the smaller grid, which has higher spatial 
resolution at the outer disc edge, in order to investigate the 
effects of tidally induced differential (retrograde) precession. 
The basic effect of differential precession on a tilted disc is to 



flatten the disc out into the plane normal to the precession 
axis (e.g. Pringle, 1992; Scheuer & Feiler 1996), and this 
seems to be the case here. Note that even the non-linear 
behaviour which occurs for large values of F± is damped by 
the imposition of forced differential precession. 

The results from the two grids do not agree on the non- 
linear development of the disc tilt for large values of . The 
large grid, extending a factor of over 200 in radius shows 
that, as was found for discs in an AGN type environment 
(Pringle 1997), the disc can become completely inverted for 
large enough values of F+. The disc then displays chaotic 
behaviour as a result of the interaction between disc self- 
shadowing, which stabilises the outer disc, and the viscous 
time-delay between the outer and inner disc radii. In con- 
trast the small grid indicates that the disc undergoes pe- 
riodic relaxation oscillations. It seems likely that this be- 
haviour is the result of too few grid points being available 
to provide an sufficiently accurate description of the disc 
shape and dynamics. The stability analysis and precession 
periods for marginally unstable discs, however, agree very 
well between the large and small grids if one accounts for 
their difference in surface density. 



4 APPLICATION TO SOME OBSERVED 
SYSTEMS 

4.1 Numerical values of the coefficients 

4-1.1 Viscosity 

As usual, the viscosity in an accretion disc is very un- 
certain. We shall make the standard assumption of an a 
disc (Shakura & Sunyaev 1973) which radiates the gener- 
ated heat locally. In the range of temperatures and den- 
sities relevant to the outer regions of accretion discs in 
bright X-ray binaries, the dominant radiative opacity is 
bound-free, k = KopT~ 7 ^ 2 . For solar abundance material, 
kq = 1.5 x 10 24 g -2 cm 5 K 7 ' 2 . Using the standard assump- 
tions we get the viscosity as a function of R and S, from 
which we eliminate E using vYl — M /3n (i.e. assuming we 
are far from the disc edge). This gives 



v = 5.; 



x 10 15 Q 4/5 M! / 8 10 M- 4 1/4 i?? 1 4 cm 2 g- 1 . 



(15) 



(M_ 8 = M/10 _8 M s yr -1 ; Mi. 4 = M/1.4M©; R n = 
R/10 11 cm). Since we measure time in units of the viscous 
time at the inner edge of the disc only the radial scaling 
enters into our calculations directly. 

When mass is added at a radius R a dd which is sig- 
nificantly less than the outer disc radius Rtide the usual 
steady-state disc assumption 3nvY, — M does not apply for 
R > -Radd- Instead, we have a regime of zero radial velocity, 
which implies vTiR^Q! =const. (Pringle 1996; prime denotes 
radial differentiation). This results in vT, oc R~ 1/2 . We can 
then compute the disc temperature by equating half the lo- 
cal dissipation rate, ^i/T,(RO,') 2 to the local emission rate 
from each disc surface. Since at its inner edge this part of 
the disc does connect to a regime where E follows from the 
accretion rate in the usual way we can still set the normal- 
isation of the outer-disc surface density from the accretion 
rate. The net result is the same as above (eq. |l^), but with 
an additional factor (_R/i? a dd) _3,/20 . We do not consider any 
discs with Rtida/ R&dd > 4, so we have ignored the small 
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correction to the viscosity in our calculations. Note that 
the surface density is significantly affected, changing from 
E oc i? -3 / 4 inside i? a dd to S a R~ 11 ^ 10 outside, so that the 
outer disc becomes significantly lighter and more suscepti- 
ble to radiative warping when mass is injected well within 



4-1.2 Disc size and mass input 

In a binary system the accretion stream emerging from the 
companion at Li initially has a non-circular orbit, but once 
it self-intersects it settles in a circular orbit with radius 
Rj which has the same specific angular momentum about 
the accreting object as the incoming stream. This radius is 
tabulated by Flannery (1975) for 0.053 < q < 19, where 
q = Mdonor/Maccrotor- We find that his calculations are well 
fitted by 



Rj 



0.085 q~ 



! (l + 0.47g 2 



(16) 



where a is the semi-major axis of the orbit. On a viscous time 
scale, the ring spreads inwards and outwards; the outward 
expansion is halted by tidal angular-momentum extraction 
when the disc approaches the Roche radius. We adopt the 
outer radii calculated by Papaloizou and Pringle (1977) as 
a function of mass ratio, which to a good approximation 
give -Rtidc = 0.87-Rl; here Rl is the volume- average Roche 
radius of the accretor, for which we use the approximation 
formula by Eggleton (1983). The inner edge of the disc is ei- 
ther near the surface of the compact object or near the mag- 
netopause, and it is at least 3 orders of magnitude smaller 
than -Rtide- Simulations with very large values of -Rtide /i?in 
are very time consuming. Setting -Rtide/ Rin = 30 — 40 in the 
numerical studies proved adequate, although we un dertake 
some calculations with Rtide/ Rin = 220 (Section 3.3). 

An important issue is where matter enters the disc, and 
here there is a clear hysteresis. If the disc is stable, mass 
enters it at the outer edge because it is stopped from fur- 
ther infall by the matter already present. But if the disc is 
tilted out of the plane, the incoming stream does not stop 
until it hits the line of nodes or reaches the point of closest 
approach for a ballistic orbit, approximately Rj/2. In prac- 
tice, therefore, there is a variation of the input radius with 
the relative phase between orbit and precession that intro- 
duces extra short-term variations in the disc. In the spirit of 
our present approximation of neglecting these rapid varia- 
tions, we add matter to the disc in a small region around Rj, 
where most of the mass would be added in a more realistic 
case. Because a flat disc with input at the edge has a signif- 
icantly higher surface density than the corresponding tilted 
disc, this means that X-ray sources can display a hysteresis 
effect: if a binary with a long period temporarily decreases in 
brightness, for whatever reason, and remains less luminous 
for a few outer-disc viscous times, it will not necessarily re- 
turn to its long period when the luminosity is restored to the 
previous level, because the now higher density means that 
a higher luminosity is needed to re-establish the warp than 
was needed to maintain the previously established one. 



4-1.3 Time scales 

There are three time scales that characterise the competing 
processes of viscous damping, forced precession, and radia- 
tive growth. In the disc regime R < J? a dd and for the above 
viscosity, their values are 



t v 



2R 2 

1*2 

R 



= 39.4 a 



-4/5 



n 



rl/4 



M~ 



tr= - = 1.13 a" 4 / 5 



e 

or 



3/10 ,,5/4 



,,3/4 ,-,-3/10 D 3/4 , 

M X ' A M_ 8 ' R^ days, 



Todays, (17) 



(18) 



2tt 



= 7.92 



lday 



1 + 5 



,1/2 D -3/2 , 



(19) 

Here e = L,/Mc 2 measures the radiative efficiency of the 
accretion process. The radiative growth time, tr, is espe- 
cially uncertain because various factors may intervene to 
make the radiation absorbed by a disc element different from 
(and usually less than) what we would estimate from our 
direct observations of the X-ray flux. In particular, the cen- 
tral source may not radiate completely isotropically. And 
furthermore, estimating the X-ray flux intercepted by the 
disc from its brightness in the UV/optical is complicated by 
the fact that only a fraction 1 — A of the incident flux is 
absorbed, where A is the disc albedo. For X rays incident 
perpendicular to a stellar atmosphere, values of A in the 
range 0.3-0.5 apply, but for the grazing incidence we have 
here it may be less. Simultaneous measurements of the X- 
ray and reprocessed optical fluxes from X-ray binaries imply 
very small fractions, /, of reprocessed radiation. For exam- 
ple, the effective temperature of the disc in CenX-4 implies 
/ = 0.07 (Heemskerk & van Paradijs 1989). In X-ray bursts 
with coincident detection of bursts of reprocessed optical 
radiation the implied value of / is at most 10 -3 (see Van 
Paradijs & McClintock 1995). It is not likely that such low 
values are entirely the result of albedo; possibly, the outer 
parts of the disc are shadowed from the central source, due 
to non-monotonic behaviour of the disc thickness. Since in 
warped discs the tilt is usually rather larger than the ex- 
pected disc thickness, this effect is probably not so impor- 
tant in our present considerations. In all, we cannot expect 
the measured precession times to agree with the simple cal- 
culations here to better than a factor few. Below we shall 
make the comparison between the observed long periods and 
expected radiative periods, where we estimate the lat ter as 
Pr = 27rfr (-Rtidc) in view of the results from Sect. 3.4. 



4.2 System parameters and simulation parameters 

In order to determine which numerical parameters are rele- 
vant to a given binary system we have to know the binary 
parameters and compute the appropriate values of various 
numerical parameters. In particular we need estimates of the 
value of 7 at Rtide, "/tide, which is a measure of the strength 
of the instability, of u> p>t idc, which is the relevant dimension- 
less measure of the forced precession rate at the outer disc 
edge, and of to and r\. In Table [l] we list the basic parameters 
of a number of X-ray binaries. Most contain neutron stars, 
for which a mass of 1.4MQ has been assumed. Parameters 
derived from the basic ones are given in Table |^. 
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Table 2. Derived parameters using the relations in this paper, assuming a = 1, 7J = 1, and e = 0.1 (and Afx = I.4M0). The time scales 
in the last three columns are taken at -Rtide- 



name 


q 






Rtide, 11 


Rj /Rtide tv 2 

(days) 


f S2„ 
p 

(days) 


Pr 

(days) 


LMC X-4 


11.4 


9.5 


0.20 


1.7 


0.26 


53. 


8.1 


7.4 


Cen X-3 


12.9 


13. 


0.19 


2.2 


0.26 


110. 


12. 


13. 


SS433 


2. 


50. 


0.32 


14. 


0.25 


890. 


100. 


110. 




14.3 


34. 


0.19 


5.4 


0.26 


800. 


47. 


62. 


LMC X-3 


0.6 


11. 


0.42 


3.9 


0.29 


250. 


21. 


62. 


SMC X-l 


12.9 


19. 


0.19 


3.3 


0.26 


110. 


22. 


11. 


Cyg X-l 


2.1 


34. 


0.32 


9,1 


0.25 


900. 


44. 


180. 


Her X-l 


1.7 


6.5 


0.34 


1.9 


0.25 


130. 


14. 


17. 


X 2127+119 


0.64 


3.1 


0.42 


1.1 


0.29 


180. 


8.7 


31. 


Cyg X-2 


0.5 


17. 


0.44 


6.6 


0.30 


420. 


140. 


29. 


gen. LMXB 


0.36 


1.2 


0.47 


0.51 


0.33 


34. 


3.3 


8.5 


X 1916-053 


0.07 


0.36 


0.60 


0.19 


0.50 


10. 


1.8 


4.0 


4U 1626-67 


0.02 


0.31 


0.68 


0.18 


0.75 


12. 


4.0 


4.9 



Table 1. Basic parameters of systems adopted in this study. 



name 




M x 


^donor 


M_ 8 






(days) 


(M a ) 


CM©) 




(days) 


LMC X-4 


1.408 


1.4 


16. 


3. 


30.4 


Cen X-3 


2.09 


1.4 


18. 


0.9 


-140. 


SS433 


13. 


10. 


20. 


10. 


164. 


X 1907+097 


8.38 


1.4 


20. 


0.05 


42. 


LMC X-3 


1.7 


10. 


6. 


3. 


198. 


SMC X-l 


3.9 


1.4 


18. 


4. 


~55. 


Cyg X-l 


5.6 


16. 


33. 


2.5 


294. 


Her X-l 


1.7 


1.4 


2.35 


0.25 


35. 


X 2127+119 


0.713 


1.4 


0.9 


0.01 


37. 


Cyg X-2 


9.844 


1.4 


0.7 


1. 


78. 


gen. LMXB 


0.2 


1.4 


0.5 


0.1 


0. 


X 1916-053 


0.035 


1.4 


0.1 


0.1 


3.8 


4U 1626-67 


0.029 


1.4 


0.03 


0.05 


0. 



Table 3. The numerical simulation parameters appropriate to 
each system in table hi using the same parameters. 



name <^ P ,tidc 7tidc rj to 



(days) 


LMC X-4 


-21. 


0.022 


7.8 


0.38 


Cen X-3 


-28. 


0.020 


7.9 


0.75 


SS433 


-27. 


0.020 


7.5 


6.3 


X 1907+097 


-54. 


0.012 


7.9 


5.7 


LMC X-3 


-37. 


0.039 


8.8 


1.8 


SMC X-l 


-16. 


0.016 


7.9 


0.81 


Cyg X-l 


-65. 


0.032 


7.5 


6.4 


Her X-l 


-30. 


0.021 


7.6 


0.94 


X 2127+119 


-66. 


0.027 


8.6 


1.3 


Cyg X-2 


-9.6 


0.011 


9.1 


3.0 


gen. LMXB 


-32. 


0.040 


9.8 


0.24 


X 1916-053 


-17. 


0.066 


15. 


0.07 


4U 1626-67 


-9.2 


0.067 


22. 


0.08 



For most of the simulations, we use a disc that extends 
from ri n = 1 to r t idc = 30. Since the inner disc does not 
matter much for our results we scale the parameters for each 
system so that Rtide in the system corresponds to the trun- 
cation radius in the numerical grid. Mass is added to the 
disc at r.j = rvideflj / Rtide (see Table ^j). The appropriate 
value of Wp.tidc follows from 

Wp,tide = £lp(Rtide)Rtide/ v l = T^p (-Rtidc)^ (-Rtidc) 



-15.63 a 



- 4/5 m- 4 1/4 m: 



3/10 



lday 



R 



11/4 



(20) 



Wp.tide depends on the details of viscosity because it mea- 
sures the precession rate, which is independent of viscosity, 
in units of the viscous time which is of course viscosity- 
dependent. In the num erical grids, Wp.tidc = ^poi'tide ■ F rom 
the simulations in Sect. 3.2.5, correcting by a factor 2 fo r th e 



higher surface densities implied by the large grids (Sect. 3.3), 
we find that tidal shear destroys the warp of a disc near 
the stability limit when u> p o is about —0.003, which means 
w P ,tide — —35. This implies that most systems have little 



enough forced precession to have their warps survive (Ta- 
ble|^), but not by much; also forced precession may some- 
what shorten retrograde long periods (Sect. 3.2.2] ). 

Similarly, 7tidc can be related to the parameters of a 
real system via 



7tido 



ti/ 2 {Rtid 

tr(Rtid 



0.0287 1- 



VM^R-^ 



(21) 



Note that 7 does not depend on details of the viscosity. 
This is because both the surface density and the viscous 
time scale are inversely proportional to the viscosity, so the 
effect that a lower-viscosity disc has a longer damping time 
for disc tilts is exactly compensated by the fact that its mass 
also increases the growth time of the radiative instability. In 
table |^ we list the values of the numerical parameters that 
follow from these relations for all systems in table [I]. We also 
list to, which is the real time that passes in each system per 
unit of dimensionless time in the simulation. 
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Figure 5. The observed long periods as a function of mass ratio. 
The solid curve indicates the expected forced precession period. 



Figure 6. Comparison of the observed long periods with the 
expected radiative precession periods. The dashed line indicates 
equality. 



4.3 Comparing forced and radiative precession 
with data 

Let us now assume that the disc radii in real binaries are ad- 
equately approximated by -Rtidc, and that the disc tilt is big 
enough that mass input occurs predominantly at Rj. Then 
we can predict what the forced and radiative precession time 
scales are for all the systems in table |§] and compare them 
with their observed long periods. 

For forced precession, a particularly simple relation en- 
sues if we divide the result by the orbital period: 



4 ^ 0.87fl L ^ 



Rorh 



-3/2 



(22) 



Since Rl/o, is only a function of mass ratio, so is the whole 
right-hand side. This presents us with a natural ordering of 
the systems by mass ratio, which is shown in Fig || The 
curve predicted by Eq. ( p^ ) is also shown, and while it is 
within an order of magnitude of the data, it does not appear 
to predict the observed values well at all. Specifically, even 
an arbitrary vertical shift of the predicted line cannot fit 
the data well, simply because no trend of the observed long 
periods with mass ratio is evident in the data, and both at 
mass ratio unity and at mass ratio 10 the observed values 
range over a factor 15. Thus a reasonable fit can be obtained 
only by omitting some of the data points (Larwood 1998). 

When we repeat the exercise to test the radiative preces- 
sion model there is not one simple system parameter with 
which the times scale, so we plot the same ordinate as in 
the previous figure versus the predicted Pr/fort>; a rather 
different picture results (Fig.^i]). There is a reasonable cor- 
relation between predicted and observed periods over the 
range of a factor 40 in Pr/Pab, with a mean ratio some- 
what below unity (note that there was no adjustment of the 
normalisation). Since there is some freedom in the choice 
of disc albedo and a parameter, we can adjust the mean 



ratio of predicted to obs erved period to be one, e.g. by set- 
ting a = 0.27 (Sect. 4.4). Then most systems come within 
a factor 2 of the curve, with the exception of Cen X-3 and 
X 1907+097, which for very similar predicted periods have a 
factor 14 difference in observed long period. (Note that they 
also have almost the same mass ratio, so forced precession 
faces the same problem of a large difference between systems 
predicted to be the same.) On the general ability to predict 
long periods, the radiative precession model therefore does 
much better than tidally forced precession. It is perhaps re- 
markable that SS 433, which is said to have a thick disc and 
for which our model should therefore not be valid, fits the 
relation rather well. One might wish to consider the possi- 
bility that the (outer) disc of SS433 subtends a large solid 
angle not because it is intrinsically thick, but because it is 
strongly warped and tilted. 

Our model can be tested further by probing directly the 
shapes of the discs using eclipse mapping and reverberation 
mapping. In Fig.|^ we show the shape of the disc as seen in 
a simulation of a fixed-shape retrogradely precessing outer 
disc that might correspond to the case of Her X-l. Note the 
flaring towards the outer edge and the near alignment with 
the orbital plane inside the mass input radius. We remind 
the reader that the range of luminosities over which these 
stationary solutions exist is quite narrow (We found it to be 
only 30%). But whether Her X-l has a constant tilt is not re- 
ally known, though large variations in the outer-disc tilt ap- 
pear not to be needed to explain the photometry. However, 
the source is known to have anomalous on- and off states or 
dips (Giacconi et al. 1973), which could well be caused by 
a variation of the disc tilt. In Fig. ^ we show a projection 
of the disc onto the sky as seen from the central source. An 
observer at infinity would be represented by a straight line 
at constant latitude as the disc precesses. When the line is 
clear of the disc, which can happen either once or twice per 
precession period, the observer sees an on state of the X-ray 
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H H 



Figure 7. The shape of a disc with mass input well inside the 
outer edge. It flares to high tilts near the outer edge and precesses 
retrogradely under the action of radiation torques, (run with small 
grid, r add = 10, r tidc = 30, F* = 0.045) 



20 - 



T3 

3 




Figure 8. A projection of the disc shape for a mildly inclined 
retrogradely precessing disc onto the sky as seen from the central 
source, (run with small grid, r add = 10, r tide = 30, F* = 0.045) 



Figure 9. Curves marking X-ray on and off states for the pre- 
cessing disc of Fig. ^ as seen for different inclinations of the binary 
orbit to the line of sight. For high inclinations only one on state 
occurs per period, but for i = 1° one observes an alteration of 
long-on and short-on states, similar perhaps to the main-high and 
short-high states of HerX-1 (Jones &: Forman 1976). (run with 
small grid, r add = 10, r tidc = 30, F* = 0.045) 




source. For low inclinations, one sees two on states of differ- 
ent length (Fig.|^), once when the line of sight passes under 
the disc (short-on) and once when it passes over it (long-on). 

By contrast, a progradely precessing disc such as per- 
haps that of CygX-2 has its highest tilt in the middle. A 
typical example is shown in Fig. When mass input in 
such a disc is not at the very outside edge the outer part 
can be going in a retrograd e dire ction at the same time, 
with a different period (Sect. 3.2.4). Under favourable incli- 



nations an observer could see both the outer and inner disc 
periodicities in the sequence of X-ray on and off states. An 
example of this is shown in Fig. jnj. 



4.4 Hercules X-l 

For the specific case of Her X-l, where a relatively stable 
cycle exists over a long time, and the retrograde precession 
has been established well from observations, we may look at 




Figure 10. The shape of a disc with mass input close to the 
outer edge. Its highest tilt occurs at the inside and the bulk of it 
precesses progradely under the action of radiation torques, (run 
with small grid, r add = 20, r tidc = 30, F* = 0.2) 



the numbers in somewhat greater detail. Using 7 cr i t = 0.1 
and Eq. (jlj]), we can establish a maximum possible radiative 
precession time scale 



27r7critii/ 2 

82.5a" 4/5 77 _1 days, 



(23) 
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star orbiting the binary has been advanced as the cause of 
the very long period (Grindlay et al. 



Figure 11. Curves marking X-ray on and off states for the pre- 
cessing disc of Fig. as seen for different inclinations of the bi- 
nary orbit to the line of sight. For high inclinations only one pe- 
riodicity is visible, but for lower ones the beat between the outer 
disc and inner disc periods is clearly visible, (run with small grid, 
r add = 20, r tidc = 30, F* = 0.2) 



4.6 Torque reversals 

Some X-ray pulsars with long and well-sampled pulse pe- 
riod records show very extended epochs in which the spin 
period decreases, followed by equally long spindown epochs. 
A famous example is 4U 1626-67 (Nelson et al. 1997, 
Chakrabarty et al. 1997), and another is GX 1+4. We note 
that in our simulations many of the more luminous discs 
are tilted through more than 90 degrees in their inner parts 
(Figure Q), implying that they are rotating counter to the 
usual direction and would provide a torque that is always 
spindown, assuming the pulsar spin is prograde relative to 
the orbit. The outer disc, at the same time, can be tilted 
much less strongly (Figure ^) . Since the range in luminosity 
over which the discs are unstable but only have small tilt 
is fairly small, and since all X-ray binaries show long-term 
flux variations, it is quite possible that a few binaries have 
luminosities that vary around the value that separates pro- 
grade and retrograde inner discs. More detailed discussion 
of this torque reversal model is given in a companion paper 
(van Kerkwijk et al. 1998). 



where we have used the parameters of HerX-1 from our 
tables. The stable long period indicates that Her X-l must be 
within 30% of the stability boundary (Sect.^J), so supposing 
we equate this period to the actually observed value of 35 d, 
this gives 



4/5 
Ol Tj 



2.36. (24) 

If we then use rj = 1/ (2a 2 ) (Kumar and Pringle 1985, Ogilvie 
1998) we find a — 0.27, which is not an unreasonable value 
(see Section 5). Note that our results imply quite generally 
that the measurement of a precession period constitutes a 
crude measurement of the disc viscosity. 

Rather than assuming near-criticality, we may also just 
compute Pr ~ 2ntr from Eq. ( [t8| ) to get 

(25) 



Pr = 17.3cr 4/5 — 



which shows good agreement with the observed 35-d period 
for reasonable values of a and e. 



4.5 Very long periods: 4U 1820 30 and 
X 1916-053 

Two well-documented long periods which are not shown in 
Figs. | and | are the 176-day period of 4U 1820-30 and 
the 199-day period of X 1916-053 (Priedhorsky & Terrell 
1984b). The reason is that both would be very far off the top 
of the vertical axis: these periods are too long in such com- 
pact binaries to be a disc precession period. For 4U 1820—30 
the bursting behaviour changes with the 176-day cycle, sug- 
gesting that indeed the accretion rate varies though the cycle 
rather than just the aspect of a disc. In X 1916—053 a third 



5 DISCUSSION 

5.1 Warp wave propagation 

The equation we use to describe the time-evolution of the 
disc tilt makes the basic assumption that angular momen- 
tum is transferred between disc elements predominantly via 
viscous processes (Pringle 1992). However, fluid discs can 
also propagate angular momentum, and therefore warps, 
through wave-like processes (Lubow & Pringle 1993, Ko- 
rycansky & Pringle 1995, Papaloizou & Lin 1995, Lubow 

6 Ogilvie 1998). For wave-like propagation to dominate, it 
is necessary that the angular velocity in the disc be closely 
Keplerian, and for the viscosity to be sufficiently small (see 
the discussion in Pringle 1997). By closeness to Keplerian 
we mean that 



|Q 



-k\ H 

«T — 

R ' 



(26) 



throughout the disc, where « is the epicyclic frequency. Mak- 
ing the approximation that the binary companion can be 
represented as a ring of mass M2 at distance a from the 
primary, mass Mi, we find that, to a first approximation in 
(R/a), ' 



5q (RV 



(27) 



Taking the system parameters for Her X-l in Tables 2 and 
3 we find that near the edge of the disc, R = -Rtide, 



|Q 



= 0.0261 



f— V 

V -Rtidc / 



By comparison we find that 



— = 0.04a 
R 



-1/10 



MZI' W MZT 



( R \ 1/8 



(28) 



(29) 
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Thus close to the disc edge we expect wave propagation 
to be marginally possible by this criterion, and to improve 
significantly at smaller radii. 

However, for wave propagation we also require that the 
disc viscosity be sufficiently small. A disc warp wave in a 
Keplerian disc propagates a distance of order H/a, before 
it is dissipated by viscosity. Thus for a wave to be able to 
propagate over a radial distance R, we require 



a < 



H 



(30) 



Measurement of a in the discs in the binary systems we are 
discussing is not easy, since the magnitude of the viscosity 
is best estimated from measuring the time scale on which 
the disc as a whole varies. However, good estimates for a 
are available for a related set of binary systems, the X-ray 
novae. For these systems it is possible to parametrise the 
quantity a by modelling the X-ray light curves through the 
outburst cycle (mainly the decline from outburst), and the 
result obtained, which it turns out also fits the outbursts of 
dwarf nova systems, is (Cannizzo, Chen & Livio, 1995) 



50 



(!) 



For the value of H/R appropriate here (eq. ES 
sponds to 



(31) 



this corre- 



ct = 0.45M" 9/46 M- 



-45/92 / R \ 
' V ^tidc / 



15/92 



(32) 



which is also in line with the values of a = 0.1 — 1.0 required 
for discs in dwarf nova outbursts (e.g. Cordova 1995), and 
for the value of a which we require to bring our calculations 
into agreement with the parameters of Her X-l (Section 4.4). 
Thus wave propagation can occur in these discs only if 



f «o.oiMr 8 1/n Mr 4 5/22 

Ft 



V fltidc ) 



5/66 



(33) 



We conclude that, in contrast to the calculations of Larwood 
(1998), the discs in binary X-ray systems are unlikely to be 
able to sustain propagating warp waves, and that the diffu- 
sion based approach is likely to provide a better description 
of the physics of warp evolution in these systems. 



5.2 Possible consequences of some neglected 
effects 

In all the simulations we used the time-independent expres- 
sion for the forced precession, i.e. we assumed the precession 
period to be much longer than the orbital period. In real- 
ity -Piong/forb is typically a few tens, but can be as low as 
5. This means deviations from the smooth precession and 
extra nodding motions of the disc with periods that are 
beats between multiples of the precession and orbital pe- 
riods can occur. These have in fact been seen in SS 433 and 
discussed in these terms (Katz et al. 1982). In addition to 
this short-period effect, mass input would also be influenced 
by the varying phase between companion star and disc: the 
accretion stream would first hit the disc somewhere along 
its curved intersection with the orbital plane, and the inter- 
section radius would vary between the outer disc edge and 
the minimum distance of approach to the centre allowed by 
the stream's angular momentum (about 0.57?j). So a more 



realistic simulation would either average the mass input over 
a large range in disc radii or have a time- variable mass input 
radius. Perhaps such effects could also explain the curious 
fact that in at least some systems (Her X-l, LMCX-4, and 
CygX-2) the precession period is an integer multiple of the 
orbital period. 

The luminosity we used in any simulation was strictly 
constant. For all stationary solutions this should be valid, 
but the solutions in which the outer disc tilt is large and 
variable (i.e. at high luminosity; Sect. 3.3) are not strictly 



valid: their surface densities are rather lower when the tilt is 
large than when it is small, and since mass is added to the 
disc at a constant rate this must mean that the amount of 
mass passing through the inner disc edge varies with time, 
and therefore the central source luminosity should do so too. 
This extra feedback operates with a delay of order the outer 
disc viscous time, and such delayed feedback could lead to 
further instabilities (as any thermostat designer knows). 

Schandl & Meyer (1994) and Schandl (1996) have stud- 
ied the effects of a disc emitting wind. Like radiation the 
wind carries momentum and therefore its launch from the 
disc surface causes a force on the disc. Since the momentum 
per unit energy in a slow wind is much greater than in ra- 
diation this could be a more important effect. However, we 
understand the formation of winds only poorly and there- 
fore the ab initio calculation of how much wind there is 
and how much momentum it carries away is difficult. Un- 
fortunately, Schandl & Meyer neglected self-shadowing and 
used an older and incorrect equation of motion for the disc, 
so it is difficult to make a detailed comparison of the two 
models and the reasons why they differ. The predicted disc 
shapes are strikingly different though, so the issue may well 
be observationally decidable: their discs have tightly wound 
spirals with the angle between the line of nodes and the 
radial direction decreasing outward, whereas the radiative 
warps have very open spirals, and the angle of the line of 
nodes with the radial direction increases outwards. If the 
sonic point of the wind were close to the disc, unlike in 
the study by Schandl & Meyer, then the behaviour would 
be almost the same as for our radiation-driven precession, 
except that the force could be much greater for a given inci- 
dent flux. This means that the region of instability moves to 
lower luminosities. Since value of 7 cr i t we find for our model 
is not too far from the actual values in X-ray binaries, this 
would mean that wind-driven warping would give extremely 
unstable, erratic behaviour in those systems, which appears 
to contradict the regular behaviour of Her X-l. 



6 CONCLUSION 

We have examined the viability of radiation-driven warping 
and precession of thin accretion discs as a model for the long 
or third periods in X-ray binaries. The model holds up very 
well and has a number of advantages over earlier propos- 
als, (i) It not only causes a disc to precess once tilted, but 
also gives rise to and maintains the tilt, whereas in previous 
models the origin of the tilt was always a difficult issue, (ii) 
The equation of motion for the disc we use does not suffer 
from the physical inconsistencies of some earlier proposed 
equations of motion. Also, unlike linearised approximations 
to the solutions of our equation in earlier papers and some 



© 0000 RAS, MNRAS 000, 000-000 



Warped accretion discs in X-ray binaries 15 



other studies, this numerical study takes account of the im- 
portant effects of self-shadowing of the disc, (iii) The quanti- 
tative agreement between the observed long periods of X-ray 
binaries and the results of our simulations is good, whereas 
forced precession due to the companion's gravitational pull 
on the disc fares rather poorly when compared quantita- 
tively with the data. Also, radiative precession can be pro- 
grade as well as retrograde, unlike forced precession. There 
is tentative evidence that the discs in CygX-2 and possibly 
X 1916—053 do precess progradely. To put it succinctly, we 
have shown that if the radiative instability of Pringle (1996) 
gives rise to the disc tilt in these systems, then it also auto- 
matically gives rise to disc precession at approximately the 
observed rate. Moreover, we have also shown that if tidally 
induced precession becomes dominant, then the instability 
is likely to be stabilised, and the disc to remain in the orbital 
plane. 

In addition to steadily precessing discs, radiative warps 
at high luminosity can also be non-stationary: their tilt an- 
gle varies periodically with time in our simulations, and in 
realistic cases with feedback between central source luminos- 
ity and accretion rate would probably exhibit non-periodic 
behaviour as well. This may be applicable to many of the 
long X-ray periods observed in nature, which are not very 
stable in amplitude and/or period. 

One particular feature of high-luminosity systems is 
that the inner disc may tilt through more than 90 degrees, 
and thus rotate counter to the normal direction. When it 
encounters the magetosphere of a neutron star it will then 
provide a strong spin-down torque, possibly explaining the 
torque reversals seen in systems such as 4U 1626—67. One 
would expect the X-ray source to be behind the disc much of 
the time when the warp is so strong, but the strongly warped 
discs have much lower surface densities, so they could be 
(partly) transparent to X rays (van Kerkwijk et al. 1998). 
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